Modeling the morphogenesis of brine channels in sea ice 
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Brine channels are formed in sea ice under certain constraints and represent a habitat of different 
microorganisms. The complex system depends on a number of various quantities as salinity, density, 
pH-value or temperature. Each quantity governs the process of brine channel formation. There 
exists a strong link between bulk salinity and the presence of brine drainage channels in growing ice 
with respect to both the horizontal and vertical planes. We develop a suitable phenomenological 
model for the formation of brine channels both referring to the Ginzburg-Landau-theory of phase 
transitions as well as to the chemical basis of morphogenesis according to Turing. It is possible 
to conclude from the critical wavenumber on the size of the structure and the critical parameters. 
The theoretically deduced transition rates have the same magnitude as the experimental values. 
The model creates channels of similar size as observed experimentally. An extension of the model 
towards channels with different sizes is possible. The microstructure of ice determines the albedo 
feedback and plays therefore an important role for large-scale global circulation models (GCMs). 

PACS numbers: 92.05.Hj, 82.40. Ck, 89.75. Kd, 47.54.-r 



I. INTRODUCTION 

Formation and decay of complex structures depend 
on changes in entropy. In the long run structures tend 
to decay since the entropy of universe leads to a maxi- 
mum and evolves into a 'dead' steady state PQ. On the 
other hand not only living cells avoid the global thermo- 
dynamic equilibrium. A. M. Turing [5] showed in his 
paper about the chemical basis of morphogenesis which 
additional conditions are necessary to develop a pattern 
or structure. For instance, cells can be formed due to an 
instability of the homogeneous equilibrium which is trig- 
gered by random disturbances. In this sense it should 
be possible, that the habitat of microorganisms in po- 
lar areas, the brine channels in sea ice, can be described 
through a Turing structure. 

The internal surface structure of ice changes dramat- 
ically when the ice cools below -23°C or warms above 
-5°C and has a crucial influence on the species compo- 
sition and distribution within sea ice [3J H]. This ob- 
servation correlates with the change of the coverage of 
organisms in brine channels between -2°C and -6°C [I]. 
Golden et AL [5] found a critical brine volume fraction 
of 5 percent, or a temperature of -5°C for salinity of 5 
parts per thousand where the ice distinguishes between 
permeable and impermeable behavior concerning energy 
and nutrient transport. According to Perovich et AL 
[6] the brine volume increases from 2 to 37 °/ 00 and the 
correlation length increases from 0.14 to 0.22 mm if the 
temperature rises from -20°C to -VC. The permeabil- 



ity varies over more than six orders of magnitude [7J. 
Whereas Golden et AL [5] used a percolation model 
we will demonstrate how the brine channel distribution 
can be modeled by a reaction-diffusion equation similar 
to the Ginzburg-Landau treatment of phase transitions. 
A molecular dynamics simulation shows the change be- 
tween the hexagonal arranged ice structure and the more 
disordered liquid water structure [8]. 

After a short introduction into the key issue of the 
structure formation we describe the brine channel struc- 
ture in sea ice and propose a phenomenological descrip- 
tion. For the interpretation of the order parameter we 
discuss some microscopic properties of water using molec- 
ular dynamics simulation in the next chapter [TTJ In chap- 
tcr |III| wc consider the phase transition and the conditions 
which allow a structure formation. We verify the model 
on the basis of measured values in chapter |IV| and give 
finally an outlook on further investigations in chapter [V] 



II. MICROSCOPIC PROPERTIES OF WATER 

A. Formation of brine channels 

Various publications report on the life condition for 
different groups of organisms in the polar areas in brine- 
filled holes, which arise under certain boundary con- 
ditions in sea ice as base- or brine channels (lacuna) 
[lOj E] . They are characterized by the simultaneous 
existence of different phases, water and ice in a saline en- 
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FIG. 1: SEM-image of a cast of brine channels [9]. 



vironment. Because already marginal temperature vari- 
ations can disturb this sensitive system, direct measure- 
ments of the salinity temperature, pH- value or ice crys- 
tal are morphologically difficult Weissenberger 
et AL [T2] developed a cast technique in order to exam- 
ine the channel structure. Freeze-drying eliminates the 
ice by sublimation, and the hardened casts illustrate the 
channels as negative pattern. Figure [T] shows a typical 
granular texture without prevalent orientation. 

Sometimes, both columnar and mixed textures occur. 
Using an imaging system Light et AL [3] found brine 
tubes, brine pockets, bubbles, drained inclusions, trans- 
parent areas, and poorly defined inclusions. Air bub- 
bles are much larger than brine pockets. Bubbles pos- 
sess a mean major axis length of some millimeters and 
brine pockets are hundred times smaller [T3]. Cox et 
AL [2J [T21 US] presented a quantitative model approach 
investigating the brine channel volume, salinity profile 
or heat expansion but without pattern formation. They 
also described the texture and genetic classification of the 
sea ice structure experimentally. A crucial factor for the 
brine channel structure formation is the spatial variabil- 
ity of salinity [17] . 

Different mechanisms are employing the mobility of 
brine channels which can be used to measure the salinity 
profile |17j . Advanced micro-scale photography has been 
developed to observe in situ the distribution of bottom 
ice algae [IS] which allows to determine the variability 
of the brine channel diameter from bottom to top of the 
ice. By mesocosm studies the hypothesis was established 
that the vertical brine stability is the crucial factor for 
ice algae growth [TIJ]. Therefore the channel formation 
during solidification and its dependence on the salinity 
is of great interest both experimentally and theoretically 
[20] . Experimentally Cottier et al [17] presented im- 
ages, which show the linkages between salinity and brine 
channel distribution in an ice sample. 




FIG. 2: Hexagonal ice (left) and liquid water at 300K (right). 

To describe different phases in sea ice dependent on 
temperature and salinity one possible approach is based 
on the reaction diffusion system 

dw(x.t) , „ , 

y-L = f(x, t) + DV 2 w(x, t), (1) 

at 

where w — ( " ~\ is the vector of reactants, x — (x, y, z) 
the three dimensional space vector, / the nonlinear re- 
action kinetics and D — f ^ ® ^ the matrix of dif- 

fusitivitics, where Di is the diffusion coefficient of water 
and D 2 is the diffusion coefficient of salt. For the one- 
and two-dimensional case we set y — z — respec- 
tively z — 0. The reaction kinetics described by f(x, t) 
can include the theory of phase transitions by Ginzburg- 
Landau for the order parameters. Referring to this M. 
Fabrizio |2TJ presented an ice-water phase transition. 
Under certain conditions, spatial patterns evolve in the 
so-called Turing space. These patterns can reproduce the 
distribution ranging from sea water with high salinity to 
sea ice with low salinity. The brine channel system exists 
below a critical temperature in a thermodynamical non- 
equilibrium. It is driven via the desalination of ice during 
the freezing process that leads to a salinity increase in the 
brine channels. The higher salt concentration in the re- 
maining liquid phase leads to a freezing point depression 
and triggers the ocean currents. 

B. Different states of water 

Already Wilhelm Conrad Rontgen had described 
the anomalous properties of water with molecules of the 
first kind, which he called ice molecules and molecules 
of the second kind |22j and which represent the liquid 
aggregate state. Dennison [23] determined the ordinary 
hexagonal ice-I modification from X-ray pattern method- 
ically verified by BRAGG [23]. This so called E h -ice is 
formed by four oxygen atoms which build a tetrahedron 
as illustrated in figure [2] In Eh-ice each oxygen atom is 
tetrahedrically coordinated by four neighbouring oxygen 
atoms, each accompagnied by a hydrogen bridge. The 
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arrangement is isomorphous to the wurtzite form of zinc 
sulphide or to the silicon atoms in the tridymite form of 
silicon dioxide. Bjerrum [55] and Eisenberg [53] have 
provided a survey about the structure differences between 
the different polymorphic forms of ice and liquid water. 

Molecular dynamics simulations with the TIP3P- 
model of water using the NAMD-software by the Theo- 
retical and Computational Biophysics Group of the uni- 
versity of Illinois show the change from a regular hexago- 
nal lattice structure to irregular bonds after the melting 
(figure [2]) . Nada et al jB] developed a better six-site 
potential model of H2O for a crystal growth of ice from 
water using molecular dynamics and Monte Carlo meth- 
ods. They computed both the free energy and an order 
parameter for the description of the water structure. Also 
Medvedev et al [127] introduced a " tetrahedricity mea- 
sure" Mr for the ordering degree of water. It is possible 
to discriminate between ice- and water molecules via a 
two-state function (G = if Mr > M T and G = 1 if 
Mr < M T ). This tetrahedricity is computed using the 
sum 



Mr 



1 



15 < I 2 > 



(2) 



1 -j 



where k are the lengths of the six edges of the tetrahe- 
dron formed by the four nearest neighbors of the consid- 
ered water molecule. For an ideal tetrahedron one has 
Mr = and the random structure yields Mr = 1. The 
tetrahedricity can be used in order to define an order 
parameter according to the Landau - de Gennes model 
for liquid crystals, which refers to the Clausius-Mosotti- 
relation. Other simulations such as the percolation mod- 
els of Stanley et al [28j [29] use a two-state model, in 
which a critical correlation length determines the phase 
transition. A mesoscopic model for the sea ice crystal 
growth is developed by Kawano and Ohashi [30] who 
used a Voronoi dynamics. 



III. REACTION - DIFFUSION MODEL 



A. 1 + 1-dimensional model equations 



We consider the reaction diffusion system 



du(x, t) 

at 

dv(x, t) 
dt 



= a\u — cu + du 5 + b\v + D\ 



d 2 u(x, t) 
dx 2 



-a-iv — b%u + D- 



d 2 v{x,t) 
dx 2 



(3) 
(4) 



in one space dimension. The order parameter according 
to the Ginzburg-Landau-theory is u(x, t) with u TO , n < 
u c < u max and proportional to the tetrahedricity u ~ 
Mt- If the variable u is smaller than u c (u < u c ), the 
phase changes from water to ice and vice versa. Thus, 
changes in u reflect temperature variations. The variable 
v is a measure of the salinity. The coefficient ai depends 
on the temperature T as (T — T c ) /T c with the critical 



temperature T c . The salt exchange between ice and water 
is realized by the gain term b±v and the loss term — 0,2V. 
The positive terms a\U and b\v are the temperature- and 
salt-concentration-dependent "driving forces" of the sys- 
tem. 
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FIG. 3: Landau-function |9| versus the dimensionless order 
parameter (tetrahedricity) for various 5. 

In order to realize the T-dependent phase transition, 
one can expand the order parameter in a power series 
corresponding to Ginzburg and Landau in equation [3] In 
order to describe properly a temperature-induced phase 
transition of second order an expression —cu 3 is neces- 
sary. The first-order phase transition is dependent on 
du° . Supercooled or superheated phases can coexist, i.e. 
an hysteresis behavior is possible. Without the term du 5 
we can also realize a brine channel formation. But this 
second-order phase transition does not allow us to con- 
sider the specific heat as a jump in the order parameter 
tp (figures [3] and [3]). We write the equations ([3| and Q in 
dimensionless form (see appendix) by setting r = \/b~b2t, 
t = f/hh/Dlx, $ = tfcyhhu, p = {Jb x c 2 lb\v and 
get 



dr 
d P {Lr) 
Ot 



M,p] 



d£ 2 



(5) 
(6) 



with ct\ = a,\j \Jb\bi, S — d\fb~\b2~/c 2 , o?2 = 0*2/ 1 \Jb\b2 and 
D = D2/D1 as well as the reaction kinetics 



f[ip,p] = axip - tp 3 + Sip 5 + p 
9bl>,p] = -a 2 p-tp, 



(7) 

(8) 



where is the dimensionless order parameter of the wa- 
ter/ice system and p the dimensionless salinity. Thus the 
dynamics only depends on four parameters a%, S and 
D. Without the salinity p in equation d3l respectively 
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^ the above equation system is reduced to a Ginzburg- 
Landau equation for the first order phase transition. 



First Order Phase Transitions 



When we neglect the salinity p, the integration of the 
kinetic function ^ yields the Landau function for the 
order parameter ip of water/ice 



F = 



2 V 



1 



4 r ' 6 

as plotted in figure [3] It possesses three minima 



(9) 



(10) 



When several different minima of equal depth exist, then 
there is a discontinuity in tp due to Maxwell construction 
and one has a first-order phase transition [31]. This is 
the case if 



and 



in consequence of 



6 = 5 r 



16a 



lc 



AS ' 



(11) 



(12) 




0. (13) 



Thus the critical parameter S = 5 C is determined by 
the temperature-dependent critical value a\ — a± c . The 
jump in figures [3] and [4] is a measure for the latent heat 
of the phase transition from water to ice. Feistel and 
Hagen [35] have deduced theoretically the latent heat of 
sea ice for various salinities. 



B. Linear Stability Analysis 

First we perform a linear stability analysis by lin- 
earizing the equation system ^ and (6j) according to 
ij) = "00 + "0 an d p = Po + P- We obtain the characteristic 
equation for the fixed points with 



4>i 



exp{\{n)T + ik£) 



(14) 



as 



¥ 



ai - 3Vo + 5 ^o 1 
-1 —oi% 



D 








(15) 




FIG. 4: The minimal order parameter ^min of ( 10 1 dependent 
on ai. 



as outlined in appendix. 

There are five fixed points for the kinetics ^ and ^ 
which satisfy / = and g — 0. In order to get a stable 
non-oscillating pattern we need a stable spiral point as 
fixed point. Moreover the associated eigenvalues have to 
possess a positive real part for a positive wavenumber, i.e. 
they have to allow to create unstable modes. Not each 
fixed point satisfies both conditions. Therefore for the 
following discussion we choose the steady state, ^0 = 
and po = 0, which corresponds to the observable brine 
channel structures measured by a casting experiment [17} 




FIG. 5: Dispersion of the linear stability (18 1 versus the 



dimensionless wave number k for a\ = 0.7, ot% = 1, D = 6 
together with the function of ( |17| l 

Short-time experiments may also record structures 
that are formed under non-steady conditions. Since those 
structures are beyond the scope of the present paper, we 
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proceed with the steady state which leads from ([15| to 

A(k) 2 + [k 2 (1 + D) + a 2 - ai]X(K) + h{n 2 ) = (16) 
with 

h(n 2 ) = Dk 4 + (a 2 - u x D)k 2 - a x a 2 + 1 (17) 
and which is readily solved 

A(k)i, 2 = ^ («i - a 2 - (1 + D)n 2 

±y/(cti +a 2 + (D- 1)k 2 ) 2 - 4) . (18) 



C. Turing space 



Let us discuss the equations (16 1 and (17 1 concerning 
the conditions for the occurrance of Turing structures in 
detail. First we concentrate on the situation k — where 
a homogeneous phase is formed. Then the fixed points 
■0 = and p = are stable according to the eigenvalues 



(181 if 



Ai, 2 (0) = - 



a 2 - a\ 



± a 2 ) 2 - 4(1 - aia 2 ) 

(19) 



are negative. Otherwise we would have a globally un- 
stable situation which we rule out. Also homogeneously 
oscillations don't describe a brine channel formation. It 



is easy to see that the solution ( 19 I gives only two nega- 
tive values if 



condition I : a 2 > ot\ and a±a 2 < 1. 



(20) 



The trajectories of salinity p and the order parameter 
tp converge to the steady state value zero by damped 
oscillations. Therefore the structure formation does not 
follow from the initial condition in time but from the 
range of the interaction in space. 

Next we discuss the spatial inhomogeneous case, k 2 > 
0, where some spatial fluctuations may be amplified and 
form macroscopic structures, i.e. the Turing structure. 
Therefore we search for such modes of ( 18 1 which grow in 



time, i.e. ReA(ft) > 0. Time oscillating structures appear 
if ImA(ft) 7^ which can be seen from the solution of ( 18 ) 
to be the case if a\ + a 2 — 2 < (1 — D)n 2 < a\ + a 2 + 2. 
In this region we have 



ReA osc (K) = 



«i — a 2 — n 2 (l + D) 



(21) 



Demanding to be positive means k 2 (1 + D) < a\ — a 2 . 
Due to (20 1 this cannot be fulfilled since the diffusion 
constant D is positive and k real. Therefore for a time- 
growing mode we do not have an imaginary part of A 
in our model. In other words we do not have oscillating 
and time-growing structures. The restriction for the only 
allowed region is 



(1 - D)k 2 -ax- a 2 \ > 2. 



(22) 



In this region we search now for the condition A > 0. The 
term before the square in (18) is negative as can be seen 



from ( 20 1 . Therefore we can only have positive A if the 



square of this term is less than the content of the root. 
This leads to 

)(a 2 + Dn 2 )>l (23) 
which restricts the k region to the interval 

^ G ( ai ° ~ a2± V(uiD + a 2 )2-4D^) . (24) 



Due to ( 20 1 the term under the square root is samllcr 



than the square of the first term in ( 24 1 and we get only 
a meaningful condition from ( 24 1 if 



a x B > a 2 . (25) 
Moreover, the square root must be real, i. e. 

( ai D + a 2 ) 2 -AD > 0, (26) 

which leads to 

n „ (1 - VI - a x a 2 ) 2 (1 + VI - U\Ol 2 ) 2 

D < s or D > s . 

af af 

(27) 



This has to be in agreement with ( 25 1 and discussing the 
different cases results finally into 



condition II 



D > 



(1 + VI -a x a 2 f 



(28) 



Having determined the ranges of a±, a 2 and D we have 
to inspect the two conditions on k, i.e. (24) and (|22| 



Discussing separately the cases D ^ 1 one sees that ( 22 1 
gives no restriction on (24 1. 



Collecting now all conditions for the occurrance of a 
Turing structure, (28), (20 1 and (24 1, we obtain 



cond.I : a 2 > ct\ and (X\a 2 < 1. 

(1 + VI - aia 2 ) 2 



cond. II 



D > 



1 



<oih1. Ill ; G [a 1 D-a 2 ± y/ (a^D + a 2 ) 2 - AD 



(29) 



The Turing space as phase diagram is determined by con- 
dition I and II and is plotted in figure [6j One can see that 
the Turing space starts at the minimal (tricritical point) 



au = ot 2t = A = 1 



(30) 



which means that we have only a Turing space for suf- 
ficient large diffusivity D > 1. For the Turing space we 
obtain the possible wave numbers according to condition 
III as plotted in figure [7] 
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FIG. 6: The Turing space as phase diagram where spatial 
structures can occur. The lower limiting line, a?2 = 
D = 1/a?, is plotted as thick line. 
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FIG. 7: The possible wave lengths k 2 where spatial structures 
can occur for D = 6 in dependence on ai and a?2- 
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FIG. 8: Dispersion h(n) for different D, ot\ — 0.7 and a?2 = 1. 

D. Critical modes 

The critical wavenumber can be found from the largest 



modes. These are given by the minimum of equation ( 17 1 
from which we find the wavenumbers 



1 



1 




and the minima 



hmin Iip9p 9^fp 



AD 



= 1 - 



(Da ± + a 2) 2 

AD ' 

(32) 
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FIG. 9: Time evolution of the order parameter ip and salinity 
p versus spatial coordinates for r = 100, 170, 400 (from above 
to below)for ai = 0.7, 0.1 = 1, 8 = , D = 6 with the 
initial condition pir = 0) = 0.5 ± 0.017V(0, 1) and periodic 
boundary conditions. 
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For n, nin > and h min 



< we find again the corre- 
sponding inequalities (|25| and (26 1. The formation of 



a spatial Turing structure, a non-oscillating pattern, re- 
quires a negative h m i n for K^in > 0. In this case there 
is a range of wavcnumbcrs which arc linearly unstable as 
seen in figure |J In figure |] we illustrate the behaviour of 
h(n) for different diffusion constants. Only those which 
lead to negative h are forming the Turing structure as 
discussed in the previous chapter. This critical range 
can be obtained, if the diffusion coefficient D is greater 
than the critical diffusion coefficient of condition II, (29) 



D r = 



(1 + VI - OtiOt 2 ) 



(33) 



which we get from h m i n = with the critical wavenumbcr 



.2 _ Defy + 9 P _ D c ai - a 2 



2D r 



2D r 



(34) 




FIG. 10: Time evolution of the order parameter ip and the 
salinity p for a\ — 0.7, a.?. = 1, S = j^- and the initial 
order parameter t/j(t = 0) = 1 and the dimensionsless salinity 
p(r = 0) = 0.5. 

The size of the structure can be estimated from ^ ZL . 
The pattern size depends on the both parameters ot\ and 
d2- The parameters determine the brine channel size 
and vice versa. With the parameters chosen in fig. [5] 
we obtain a pattern size of 12.6. In the next chapter we 
compare this value with experimental quantities. With 
a small initial random perturbation we plot snapshots 
of the time evolution of the order parameter tp and the 
salinity p in figure [9] The quantities ip and p are opposite 
to each other; domains with low salinity correspond to 
domains with ice and domains with high salinity to water 
domains. We see the formation of a mean mode given by 
the wave length n c . 



A positive h(n — 0), respectivly a negative A(k = 0), 
for k — guarantees that p and ip converge to the sta- 
ble fixed point p = 0, ip = 0. Therefore the structure 
formation does not follow from the initial oscillations in 
time (fig. 10 1. In order to obtain a new spatial structure 



there must exist at least a negative h for n > 0, respec- 
tively a positive eigenvalue A, for k > as was discussed 
in the last chapter. 



E. 2+1-dimensional model 

From the characteristic equation in the spatially two- 
dimensional case 

A 2 + + k 2 )(1 + D) + a 2 - ai}\ + h{n\, fc 2 ) = (35) 

we find the corresponding dispersion relation 

/i(k 2 , K 2 ,) = D(K^ + K.^) 2 + (a2-aiD)(K^ + K^)-aia2 + l 

(36) 

which is illustrated in figure [TT] 

The Turing space is bounded by the sectional plane 
h = 0. The evolution of the order parameter ip and 
the salinity p is illustrated in figures [12] and [13] Their 
behavior is inversely proportional and corresponds to the 
fact, that a high salinity occurs in the water phase and 
a low salinity in the ice phase. Similiar as in the one- 
dimensional case we see the dominant formaton of one 
wavelength. The model kinetics generate brine channels 
of similar size. In order to obtain a hierarchical net of 
brine channels of different size, the kinetics in the basic 
equations can be altered accordingly |33J. 
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FIG. 11: Dispersion of the two-dimensional characteristics 
(35) and ([Ml for ai = 0.7, a 2 = 1, D = 6. 
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Note on second and thirth order kinetics 



If we replace the kinetics ^ by 

f[ip,p] = a!ip-ip 3 + p 

fbP,p) = a 1 if;-ip 2 + p 



(37) 



(38) 



we can carry out the same linear stability analysis for the 
fixed point ipo = and po = 0. Then we obtain the same 
characteristic equation as (16) and (17). Therfore we 



get the same Turing space for the structure formation. 
Consequently, both kinetics allow us to realize a brine 
channel formation but the thirth order kinetics describes 
a second order phase transition only. In this connection 
it is possible to discuss second order phase transitions 
with spin models, too. 



IV. CONNECTION TO EXPERIMENTAL DATA 

The critical domain size is determined by the equation 
( 34 1 . Due to this relation we can infer other parame- 
ters in the model equations ^ and Q. From the re- 



lation between the dimensionless wavenumber n and the 
dimensional wavenumber k 



2D 



Pi 

Vhb- 



:k 2 



we get 



2tt 



12.6 



2tt Vhh 



(39) 



(40) 



The observed diameters of the brine channels range 
from pm to mm scale [9]. For a size of 2tt jk c — 10 pm 
and a diffusion coefficient D\ = 10 -5 cm 2 s -1 for H2O- 
molecules we obtain the product &1&2 = 2.5 • 10 6 s -2 and 
a transition rate a\ = \fb1b2a1 = 1111s -1 . The rate ai 
is proportional to reorientations of the molecules per sec- 
ond, l/r d = 10 5 s -1 (Eisenberg) [53] and to the scaled 
temperature 



a 1 



T c -T 1 



(41) 



where T c is the melting point depending on the salin- 
ity. The mean salinity in sea ice of 35g/l corresponds 
to 1 7VaC7-molecule per 100 ^O-molecules, i.e. 1 
Na + -ion and 1 C7 - -ion per 100 ^O-molecules in a di- 
luted solution after the dissociation or a ratio of x = 
{ n Na+ + n ci-)/ n H 2 o = 1/50. From this facts we obtain 
according to Clausius-Clapeyron 



AT 



xRT 2 
AH 



(42) 



a freezing point depression from 0°C to — 2°C, where 



AH = 6-^t is the latent heat of the phase transi- 



tion from water to ice, R 



.314 



is the univer- 



rnolK 

Thus we obtain cor- 



sal gas constant and T = 21ZK. 
rectly T c = 271K . For an environmental temperature of 
T 



-5°C 



268-R" according to (41 1, a transition rate 

s~ 

1111s -1 , which we estimated from the 



of 271-268 • 10 5 s -1 = 1107s -1 follows which nearly cor- 
responds tO CL\ 



domain size ( 39 1 
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FIG. 12: Structure formation for 3 time steps r = 
100, 170, 400 (from top to bottom) for the order parame- 
ter ^ (left) and the salinity p (right). The parameters are 
ai = 0.7, 02 = 1, d = ttt — , D = 6 with the initial con- 
dition v(t = 0) = 0.5 ± 0.01iV(0, 1) and periodic boundary 
conditions. 
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FIG. 13: Magnification of a detail of figure [T2"| 

Furthermore we find the transition rate a 2 = 
V / 5i&2«2 = 1587s -1 and the diffusion coefficient D2 = 
6 x 10 -5£IS -. Due to the transformation of equation &3\ 
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into the dimensionless form ([5| there exists a fixed rela- 

We obtain for the equation 



tion (11) because of c = 1 . 

([3]) the rate d c with a\ = a,\j \Jbib 2 , S = d\Jb\b 2 jc? of 



d c 



3c 2 
16ai ' 



(43) 



which is proportional to c 2 . A transition rate c = 
1000s -1 yields a critical rate d c — 169s" 1 . From the 
knowledge of the diffusion coefficient Z?i and the size of 
brine channels we can deduce the two rates a± and a 2 . 
Both rates possess the same order of magnitude and are 
inside the Touring space of structure formation. If the 
experiments would lead to other parameters a\ and a 2 , 
especially to rather different values, a brine channel could 
not arise because of the limitation of the Turing space in 
figure [6] and [7j In other words the model here seems to 
describe the experimental finding of brine channel forma- 
tion. 

Due to the small difference between the time constants 
a\ and a 2 we obtain a dynamic interference between the 
reorientation of the water molecules and the desaliniza- 
tion. Both are evolving on nearly the same time scale. In 
particular, we cannot simplify the kinetics by separating 
time scales using the Tichonov theorem [34 in order to 
reduce our reaction diffusion system (J3J , ([4j) but have to 
consider both dynamics as demonstrated. 



V. CONCLUSION 

In this paper it has been shown that a reaction dif- 
fusion system which connects the basic ideas both of 
Ginzburg and of Turing can describe the formation of 
brine channels with realistic parameters. For the cho- 
sen parameters patterns of similar size emerged. Eicken 
[35] and Weissenberger [9] distinguished between six 
various texture classes of sea ice dependent on the crys- 
tal morphology, brine inclusions and the genesis. The 
different ice crystal growth depends on snow depostion, 
flooding, turbulent mixing, quiescent growth rate or su- 
percooling. Each condition determines the character of 
the kinetics. Non-linear heat and salt dissipation for ex- 
ample leads to dendritic growth (snowflakes) whereas one 
observes in sea ice mostly lamellar or cellular structures 
rather than complete dendrits [7J. Hence, the morphol- 
ogy of sea ice is one criterion for the choice of an appro- 
priate kinetics for the genesis of sea ice. Therefore, in 
order to simulate different structure sizes and textures, 
we can modify the dispersion relation by varying the pa- 
rameters ax, a 2 and D or by a modified kinetics |36j . The 
crucial point is the shape of the dispersion function. If 
there are multiple different positive unstable regions for 
the wavenumbers with positive real part of eigenvalues we 
could expect that differently large channels evolve. For 
instance Worster et AL [20] has presented a general 
theory for convection with mushy layers. The two differ- 
ent minima of the neutral curve, determined by the lin- 
ear stability analysis, correspond to two different modes 



of convection, which affect the kinetics and determines 
the size distribution of the brine channels. We note that 
the initial conditions are decisive for the appearance of 
specific pattern [33]. Hence one should investigate how 
dislocations or antifreeze proteins influence the formation 
of the brine channel distribution. 



APPENDIX A: DIMENSIONLESS QUANTITIES 



If we set t — -r- , £ 



with 



and 



to 

dip dr dtp 



at 



dt dr 



x 

if* and 

t OT 



u = C\ip and v = C 2 p we get 

di dip _ j_a^ 



dx <9£ xq d£ 



Oil 

dt 



dip Ci dip 



dt 



to dr 



du dip C\ dip 
dx 1 dx xq dt; 



Because of 2J- = §4 
f^r = ^2 fg£ and consequently 



dip d 2 x 
dx d£ 2 



r 2d 2 ip 
L dx 2 



(Al) 



(A2) 



we obtain 



d 2 u _ c d 2 ip _ Ci d 2 ip 
dx 2 1 dx 2 x 2 9£ 2 



(A3) 



Accordingly one has 



dv 
di 
d 2 v 
dx 2 



Chdp 
t dr 

x 2 di 2 ' 



(A4) 



From ([3]) and (EJ follow the dimensionless equations 



dip , r , h d 2 ip 
- = f[ip >p] + D 1 - lw 



dp 



n 1 _l n *° 8,2 P 
= g[iP, p }+D 2 ^ w 



(A5) 
(A6) 



with 



C 2 



cutoip - cf CfV> + dtoC^ + hto-^p 

<~>1 



Ci 



9bP,p\ = -a 2 t p-b 2 to—ip. 

<^2 



(A7) 



If we choose 



we obtain C\ 




L-2 



(A8) 



and 
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APPENDIX B: LINEAR STABILITY ANALYSIS 

Let ip and p denote small displacements from the equili- 
brium values ipo and po and write ip = tpo + ip and p — 
po + p- With respect to ([5j |6l we obtain 



+ S(^o + $) 5 + Po + P + ^d^ T ) 
Pt(£,t) = Dp 66 (Z,T)-a 2 {p + p)-{iPo+ip)(Bl) 

If we consider only linear terms 

$t(€,t) = h - 

Pr(^r) = ■■■-i)-a 2 p+--- + Dp ii {^T) (B2) 



we get 



— ^PQ,P — PO) 



9£ 2 
D 





« 2 



where J ( 



is the Jacobian 

_ ( ftp fp 




(B3) 



(ip=ipo,P=Po) 
J(i/>=ipo,P=Po) \ ,, I 

ai - 3iP 2 + 55^o 1 

-1 — C*2 

which we can calculate also considering ([7| 



(B4) 



Using the Fourier ansatz ip = ipoexp(\r + in^) and 
p = poexp(Xr + iK^) in (15) we find 



ay - 3^ 2 + i \ / ^ 

-1 -Ot2 ) \ Pi 

-K 2 \f & 

-Dk 2 J I pi 



With ^o = and po = the eigenvalue equation 



= 



n 2 — ai 



-1 



A 
A 



1 a 2 + £>K 2 
follows with the characteristic equation 

1 







OL\ — K 2 — A 

-1 -OL2- Dk 2 - A 
A 2 + [K 2 (l + D) + a 2 -ai]A 
I?k 4 + (eta — a\D)n 2 — a\ct2 + 1 

/l(K 2 ) 
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